clear; close all; addpath('../utilities');

para = [0.046865000000000  10.281160000000000   0.038592000000000  20.693224000000001   3.594894000000000   0.112565000000000   0.146020000000000   7.957018000000000];

main_para

options = optimoptions('fsolve','Display','iter','TolFun',1e-12);
beta    = fsolve(@solve_beta,0.99,options,gamma,psi,mu_C,sig_C,P,P_nodis);

main_SDF

lamS    = ones(nS,nK);
lamD    = ones(nS,nK);
k       = linspace(k_min,k_max,nK)';
k_bar   = ones(nS,1)*5;
k_iss   = ones(nS,1)*10;

[kgrid,mgrid]   = meshgrid(k,(1:nS)');
kgrid           = kgrid(:);
mgrid           = mgrid(:);

Q_big   = repmat(Q,nK,1);
Rf1_mat = repmat(Rf(:,1),1,nK);

main_VFI

main_CDS

main_CDS_P

money       = ones(nS,nK); % moneyness levels

main_options

IV          = NaN(nS,nK,2);
IV(:,:,1)   = BSIV;

money       = exp(-BSIV/sqrt(12)); % moneyness levels

main_options

IV(:,:,2)   = BSIV;
IV_ATM      = IV(:,:,1)';
IV_ATM      = permute(IV_ATM,[2,1]);
IV_Skew     = IV(:,:,2)'-IV(:,:,1)';
IV_Skew     = permute(IV_Skew,[2,1]);

main_Ret

main_simu

save model_simulation table_10 G_m
